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Abstract Non-linearities play an important role in micro- and nano- electromechan- 
ical system (MEMS and NEMS) design. In common electrostatic and magnetic ac- 
tuators, the forces and voltages can depend in a non-linear way on position, charge, 
current and magnetic flux. Mechanical spring structures can cause additional non- 
linearities via material, geometrical and contact effects. For the design and operation 
of non-linear MEMS devices it is essential to be able to model and simulate such 
non-linearities. However, when there are many degrees of freedom, it becomes dif- 
ficult to find all equilibrium solutions of the non-linear equations and to determine 
their stability. In these cases path following methods can be a powerful mathematical 
tool. In this paper we will show how path following methods can be used to deter- 
mine the equilibria and stability of electromechanical devices. Based on the energy, 
work and the Hamiltonian of electromechanical systems (section [U, the equations 
of motion (section |2]i, the equilibrium (section (3]) and stability conditions (section 
|6]l are derived. Examples of path following simulations (section|4|i in Mathematica 
(section|5]l, Matcont (section|7]i and using FEM methods in COMSOL (section[8ll 
are given to illustrate the methods. 
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1 Energy and work in electromechanical systems 



MEMS and NEMS devices can often be treated using classical physics. Interac- 
tions between the electric and magnetic fields can usually be neglected in MEMS 
and NEMS devices, because the system dimensions are small compared to the elec- 
tromagnetic wavelength and are large enough to neglect quantum effects. In these 
approximations the total internal energy of the system Utot and total work 

done by external sources Wtot = can be expressed as a sum over different types 
of energy and work, given by equations (fTlfTOli: 



u„n = I \^dr = W— (1) 

Jr 2 p,„ 2 y mi 

U,n,v = i VGp.dr = ^ ^ (^mt £ ^m^ (2) 

ij.Adr=ii:(^'i:A7'^.) (3) 

[,\vp,dr=^-l^(^Q,Y,Cr/Q^ (4) 



Ustrain = j ^6 ' OdY = \^iY,kijx/j (5) 

Win = / / ^ext- dp,„dr = V / Vext.idpm,i (6) 

Wgrav = j^^ VG,e,,dp^dr = T. Jj mgext.idxi (7) 

Wn,ag = [ f^iext ' dAd = £ f^' hxt,A^i (8) 



W,, = j^' VextdpedY - j^' VextjdQi (9) 

Wstrain = [ f O ex, ■ dgd^ = £ T' Fext.idXi ( 1 0) 

These expressions are given to show that there is a striking similarity between 
energy and work equations in different physical domains. As a result of this, their 
physics can be analyzed by the same mathematical methods. Moreover, physical 
coupling and transduction between domains, which are often encountered in elec- 
tromechanical systems, can also be analyzed. In this paper these mathematical meth- 
ods, based on the conservation of energy {Uu,i ~ W,oi), will be discussed and some 
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examples will be given of how numerical path following can be used to analyze such 
systems. 

The energy and work expressions in each of the equations can be written in terms 
of volume 'f integrals over energy densities, which can be treated numerically by 
finite element methods. The work equations also include an inner integral, which 
expresses the fact that the work done can be path dependent. If the system consists 
of discrete elements, simplified energy expressions on the right side of ([TlfToli can 
be given, which can substantially reduce the number of degrees of freedom and 
are more suitable to treat systems of particles or electrical networks with lumped 
element components. 

The variables in equations ([TlfTOli are given as follows. Pm is the mass density, p,„ 
is the mass momentum density. Yq is the gravitational potential and G is the grav- 
itational constant. Y is the electric potential, pe is the charge density, J is the total 
current density and A is the magnetic vector potential, a and e are the mechanical 
stress and strain. The expressions on the right side of the equations are written as 
sums over discrete charges Q,-, magnetic fluxes and masses m,- at positions x,- 
and with mass momentum p,,, which are coupled via effective springs kij, grav- 
itational potential G/rtj, capacitors C,, and inductors L,y. The components of the 
inverse capacitance and inductance matrices C^' and L^' relate potential to charge 
Vj = Y^iCfj^Qi and relate current to magnetic flux Ij = Y.iLJj^<Pi. The total work 
Wtot supplied to the system is a sum of the work from external gravitational Vc-eir' 
current density Jext, electric Vext, magnetic A^^t and stress Gext fields. The kinetic 
work term Wjt„, is needed to account for a relative speed difference \ext between the 
frame of reference of the system and the frame of the observer The work can also 
be performed on discrete elements by external gravitational acceleration gext.i, cur- 
rent Iext,i or voltage Vext.i sources or external (stress) forces Fgxtj- For more details 
on equations ([TlfTOli the reader can refer to standard texts on classical mechanics and 
electromagnetics. 

External thermal sources at temperature T^xr, which add heat to the system, 
can be treated by similar methods, although this requires keeping track of the en- 
tropy S and temperature T in the system. According to the first law of thermody- 
namics {dUtot = SQij + 5Wtoi), this requires adding an integral Qi, = Jq T^xtdS to the 
work equations. The internal thermal energy of the material is in essence a sum of 
the kinetic energy of the atoms which can be expressed by a thermal energy term 
Utherm,kin = c{T)p,„dTdy , where c{T) is the specific heat capacity. 

The work provided by an electrical source in equation (|9]l can be rewritten us- 
ing integration by parts: Wei = jr {VextPe — Sq"' PedVext)d'V . It is often convenient 
to define the coenergy or complementary energy as the Legendre transform of the 
work such that t/* = VextPedf - Wei = U lo'" PedVextd^. Figure [U shows how 
energies and coenergies can be represented graphically for a voltage controlled ca- 
pacitive MEMS switch example which is treated in section|5]and figure|3] 
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Fig. 1 Energy and coenergy in a voltage controlled capacitive MEMS switch. The areas between 
the lines in this graph represent the internal electric energy U^i and strain energy Ustmin (between 
blue dashed and black solid lines), the kinetic energy Uki„ (between red dashed arrows and black 
solid lines) and the path-dependent electric coenergy (7*, of a capacitive MEMS switch. The black 
curve will be calculated in section |5] and represents the static equilibria of a capacitive MEMS 
switch which is actuated by a voltage V„/ . If the device follows the black path, no kinetic energy is 
generated because the device is always in static equilibrium. The red dashed arrows represent the 
hysteresis path that is followed when the switch is actuated using voltage control. In this case the 
areas between the red and black lines is the generated kinetic energy during pull-in (Ukin.duse) and 
pull-out (Uki„.open)- The sum of both energies Uiti„^close + Uiti„j,pe„ is the dissipated energy during 
one switching cycle. 



2 Equations of motion and energy conservation 

The Hamiltonian is the difference between the total internal energy Utot and the 
total supplied work Wu,t ■ It can be described by a set of generalized coordinates qi 
with (/ = 1 . . .A^) and corresponding generalized momenta p,: 

(g'l, . . .,qN,pu- ■ ■ ,Pn) = Utot - Wtot (11) 
From the Hamilton equations the generalized forces Fi and speeds v, are given by: 

Fi = Pi = , Vi = qi = — — (12) 

dqi dpi 
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The total internal energy in the system equals the total work W,o, done on the 
system and therefore the Hamiltonian Jif , which is the difference between internal 
energy and work = JJtot — Wtot, is conserved. This follows from the Hamilton 
equations (fT2] i because: 

Let us choose as generalized mechanical coordinates the position qmech.i = 
with as corresponding generalized momentum p,„ech.i = Pm.i, and for the generalized 
electromagnetic coordinates the charge qem.i = Qi with the magnetic flux pe,„ j = 
as corresponding generalized momentum. 

If the parameters kij, LJj ' and Cfj ' are constants, applying the Hamilton equations 
(fT2b to equations (fTlfTTI) gives the equations of motion: 

Fmech.i = Pin.i = ^Y^kijl j + niigext.i + F^xt^i , Vmech.i^Xi = p,„,i/mi (14) 
J 

Fem,i = <P, = -Y,Ci-/Qj + Ve„J,Ven,.i = Qi = Y.^7/'^J+tex,,i (15) 

j J 

This gives us Newton's second law, Hooke's law and Faraday's law of electromag- 
netic induction. 

However, if for example the spring constants, inductance and capacitance pa- 
rameters kij, LJj^ and Cjj^ depend on the mechanical positions x,-, which is often the 
case in electromechanical systems, the mechanical and electromagnetical problems 
become coupled and the mechanical forces become: 

J l.J ' IJ ' 'J ' 

+migext.i + Fext,i (16) 

This shows that if the parameters depend on x, the generalized mechanical force 
Fmech.i becomes the sum of coupling forces from different physical domains, like 
the spring force, electrostatic and magnetostatic forces. In equilibrium this sum of 
coupling forces is zero. 



3 Electromechanical equilibrium 

The system is defined to be in electromechanical equilibrium if all generalized 
forces Fi are zero. In equilibrium, the gradient of the Hamiltonian is therefore also 
zero as shown by equations ( |17|18t : 
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Fi = pi 

F = -Vq^^ 

For a specific set of generalized momenta p, and generalized external parameters 
^gen = {n'tgextiFextiyextilext}^ equation ( fTTl ) gives generalized force equations 
with unknowns which can be solved to obtain the equilibrium state of the system. 
If the energy function has only one variable parameter and all other parameters 
are fixed, the equilibrium state ({eq.Q at ^ = thus satisfies; 

F(q.<,,o,^o) = -Vqjr(q,^,o,^o) =0 (19) 

Often equilibria exist for which all generalized speeds v,- are zero. These equi- 
libria are called static equilibria or stationary states. The equilibrium solutions of a 
non-linear system can be found using a powerful mathematical technique, which is 
called path following or numerical continuation HHH. 



(17) 
(18) 



4 Path following and numerical continuation methods 

Finding all equilibrium states of an electromechanical system is a difficult prob- 
lem that cannot be solved in general because it requires one to check the whole 
dimensional coordinate space for solutions of equation WT\ . for each value of the 
external parameters ,'^j,en and momenta p,-. However, when an equilibrium solution 
of the equation ( [TtI i is known for a certain set of parameters p,-, J^t.^,,, it is possible 
to efficiently find additional solutions using numerical path following techniques 

ma. 

Recently it has been shown that path following techniques can be very useful 
for the analysis of electrostatic MEMS devices with non-linear contact forces ifTOl . 
The technique can in fact be applied to almost any non-linear electromechanical 
system described by equations ([TTfTOli. In MEMS and NEMS devices, the system is 
usually controlled by a single external parameter which can for example be the 
voltage or current from an external source for electromechanical actuators, but can 
also be an acceleration or gas pressure force in sensor systems. The 1 -dimensional 
external parameter ^ with the A^-dimensional coordinate space q form a (A^ + 1)- 
dimensional space. 

The implicit function theorem implies lHHD that i/one equilibrium solution q^^ o 
with F(qg£y.o, ^o) = exists at external parameter value =^0j ^hen it has to be part of 
a continuous 1 -dimensional equilibrium solution curve qe<;(^) at parameter values 
^ in (A^+ 1) -dimensional space with q^^ o = qc</(^o)- It is therefore often possible 
to start at a known equilibrium solution qeq{^o) and determine the other solutions 
by following the curve iqeqi'^),'^)- 
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Path following methods usually employ a predictor-corrector algorithm HI to 
follow an equilibrium curve. In the predictor step, a new solution is predicted, usu- 
ally by taking a step in direction of the tangent of the solution curve. In the cor- 
rector step, the predicted point is brought back on the solution curve, usually by 
an iterative method. For the predictor step it is necessary to determine the tan- 
gent direction to the solution curve. Let us assume that there exists an infinites- 
imal vector tangent to the solution curve dq£.^.o(d^o), with a corresponding in- 
finitesimal parameter change d^o such that there is another equilibrium state at 
{qeq{-J^),^) ~ (q«/,o +dq«/,o,=^o +d,J^o)- Because the new state should also sat- 
isfy equation ( fTSl l the predictor algorithm can determine the tangent direction vector 
(dqe<^,0, d^o) in + 1 dimensional space by solving the equations: 

F(qe^,o + dq,^,o,^o + dJ^o) - F(q,,,,o,^o) 

= h{q,,,o,^Q)^<leq.O -H^(q^^_g_,^g)dq,^,o = (20) 

Where JF(q„,o.-^o) Jacobian of F at (qety.o,^o)- HjT is the Hessian matrix of 
the Hamiltonian, which is the Jacobian of the gradient of the Hamiltonian. 

Besides the tangent direction, the predictor algorithm also needs to specify the 
steplength of the vector. The curve {qeii{s),^{s)) can be parametrized by a param- 
eter s, such that the steplength is given by \{dqec,/ds,dJ^/ds)\ x As. If parameter 
continuation is used, the parameter s is proportional to the external parameter such 
that s = aJ^, this however leads to problems at folds for which the steplength be- 
comes infinite because Idq^g/ = °°. In the commonly used pseudo-arclength 
continuation algorithm, the parameter s is parametrized by the length along the curve 
s = Jq a\{dqgg/ds,d,^/ds)\ds such that the steplength is always finite. 

Note that in order to use path following methods to obtain equilibrium solutions, 
it is required to start from an initially known equilibrium solution. It is therefore 
essential that an initial solution can be found, this can sometimes be difficult since 
there is no general method to find these initial solutions. 



5 Example of path following method 

As an example of the methods discussed in sections [T]|4] we consider the problem of 
finding the equilibra of a voltage controlled capacitor connected to a spring from its 
energy and work equations. A schematic picture of this system is shown in figure|2] 

The system has two generalized coordinates x and Q. The internal energy and 
work of the system can be used to derive the Hamiltonian: 

^ = U,,rain + U,l - W,l = ^kx^ + ^jf— ^ " 1^ Vext .eq(.Q)dQ (21) 

2 2Aeo Jo 
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Q 




Fig. 2 An electrostatically 
actuated capacitor of which 
one plate is connected to a 
spring. 



C=Aso/(g-x) 




Note that for the conservation of energy to be vaHd, it is necessary that the externally 
applied voltage Vext,ei]{Q) is taken such that the system is always in equilibrium 
along the path, otherwise the kinetic energy of the system would become non-zero. 
The generalized force vector F is found from the gradient of the Hamiltonian: 

F = -V,^ = - — — , — — 
V dx dQ J 

The Hessian matrix of the Hamiltonian is founcQ to be: 




- = Jf = - 1^ ^i"^ J I = I o dv„,_,„ ° I (23) 

Thus equation (|20] | gives us: 

-H^(.x,.e) = (-M.+|-.e,^..+^.e-^^e) =0(24) 

Dividing these two equations by dQ and rearranging shows that: 

dj Q 



dQ kAeo 

dVexr,eq _ dVexr.eq _ 1 Q^ 



dQ dQ C kA^-el 



(25) 
(26) 



At an equilibrium point at coordinate {xeq,Qeq,yext,eq) equations dZSl l and ( |26] | 
provide the tangent vector to the solution curve {^dQ^dQ, —f^^dQ). In figure[3] 



' The method presented here is slightly different from the method which we presented in |10|. In 
1101 the electrical generalized force equation Vexi =Q/C was eliminated before taking the Jacobian 
of the force vector. In this derivation we keep all generalized force equations to have a square 
Hessian matrix. 
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Fig. 3 This script in Mathe- 
matica shows how the — Q 
curve of an electrostatically 
actuated MEMS switch with a 
non-linear exponential spring 
contact can be obtained us- 
ing path following as de- 
scribed in section |5] The 
spring energy is given by 
jkx^ + exp (ci (x - g + 5)), 
the electrical energy is 
jlC and the capaci- 
tance C = Aeo/(g — x), with 
fe = g = Aeo = \,c\ = 100 and 
5 = 0.1. The path is followed 
starting from the solution 
X = V = e = 0. The V„, - Q 
curve is identical to that in 
figure [T] The x — Q curve is 
also plotted. 

we show how the equations in this section can be implemented in a Mathematica ||8] 
script to determine the Y ~Q curve of an electrostatically actuated MEMS switch 
with a non-linear exponential spring contact. In this example we have used the ND- 
Solve function in Mathematica to follow the solution path solely from the tangent 
vector function, without using a corrector method. The different types of energy in 
this system are shown in figure [T] 

Dedicated interactive numerical path following packages Q like Matcont Q 
and Auto (|2l include coiTector methods and stability analysis. They are therefore 
more robust and convenient to treat complex problems. An example in MATCONT 
is discussed in section|7] 



6 Stability 

The system is in a static equilibrium state f\eq when the sum of all forces is zero 
and all generalized momenta are zero. This equilibrium is only stable if for all pos- 
sible infinitesimal displacement vectors dq^^, the Hamiltonian increases such that 
+dq£.(^) > Jif{qeq)- Any displacement would therefore violate the conser- 
vation of energy and thus the system will remain in state q^jg forever 
The local Taylor expansion of the Hamiltonian is given by: 



Htot = 1/ 2x^2* Exp [100 (X- 0.9) ] +Q^2 (1 -x) / 2 
-Integrate [V[Q1] , {Ql, 0, Q)] 

1 

gioo (-o.9«) ^ _ q2 _ jj, ^ V[Q1] dQl 

2 2 Jci 

Hessiaranatrix = -D[Htot, {{x, Q}, 2}] / . x -> x[Q] 

{ {-1 - 10 000 1-".S«"Q1> _ q]_ ;q, -1 + x[Q] + V [Q] ) } 
Dotprod = Hessiaranatrix .{x'[Q], 1} 

{q + (-1 - 10 000 e-°° l-0-S™l01i ) x'[Q] , -1 +x[Q] +V'[Q] +Qx'[Q] } 

NumericalPath = HDSolve [ { Dotprod [ [ 1 ] ] == , Dotprod [ [ 2 ] ] == , 
V[0] == 0, x[0] == 0), {V[Q], x[Q]), {Q, 0, 10)]; 

Plot[(V[Q] , x[Q]) /. NumericalPath, {Q, 0, 5)] 




1 2 3 4 5 



(q + dq) = (q) + VJT (q) • dq + ^dq • H^(q)dq (27) 
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At a static equilibrium point, the gradient of .yP is zero according to equation 
( [19] ). so the system is stable if the right term of equation ( |27] i is positive for all dq. 
The vector dq can be written as dq = Lf^flA/./q//,;, where q// , are the eigenvectors 
of the Hessian matrix with corresponding eigenvalues A,-. Therefore in equation dZTl i 
the term jdq • Hjf (q)dq = \ ^i'^dq All r always positive if all eigenvalues A,,- 
of the Hessian matrix are positive. Points at which one eigenvalue becomes zero 
are called folds. Points at which more than one eigenvalue become zero are called 
bifurcations. Since these zero crossings cause a sign change of the eigenvalues, the 
stability of the system usually changes at these points. 




Fig. 4 Left: Model equations in Matcont for two identical switches with the same parameters 
as the switch in figure [3] Right: Calculated equilibrium curves. Folds (LP) and bifurcations (H) 
are automatically detected by the software. The folds (LP) correspond to the pull-in and release 
instabilities of one switch, at bifurcations (H) both switches are unstable for pull-in or release. 



7 Numerical path following in Matcont 

As an example of the use of MATCONT for path following, we take two identical 
uncoupled switches with the same parameters as in figure |3] and enter the force and 
momentum equations ( fT2] i in the system window shown on the left side of figure 
|4] To improve the convergence of the calculation a mass and damping term have 
been added to the force equations. On the right side of the figure the calculated path 
following curves are shown. Path following programs like MATCONT monitor the 
eigenvalues of the Hessian on the equilibrium path to determine the stability of the 
system and plots the folds (LP) and bifurcation points (H) on the graph (see figure 
|4]i. As an example two identical capacitive MEMS switches are simulated in in 
figure|4] All eigenvalues of the Hessian are positive on the line from (0,0)-H. At the 
bifurcation point H (5,5), two eigenvalues become zero and both switches become 
unstable towards pull-in. On the paths H-LP only one switch is unstable and the 
other is open (one negative eigenvalue), on the path H-H both switches are unstable 
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(two negative eigenvalues). On the segments LP-LP between the release and pull-in 
point of one of the switches and on the segment H-(l,l), where both switches are 
closed, the system is stable and all eigenvalues are positive. 



8 Finite element method path following simulations 

The examples given up to now have a low number of degrees of freedom, such that 
their energy equations can be entered manually. To model continuous systems, for 
which the energy and work equations are given by the volume integrals in equa- 
tions ([TTfTOb. dedicated finite element method (FEM) software is available. These 
packages mesh the volume and construct partial differential equations that provide 
the equations of motion for all degrees of freedom. Moreover they provide efficient 
solvers, such that systems with thousands or even millions of degrees of freedom 
can be solved. 
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Fig. 5 Simulation of the equilibria of a voltage actuated single clamped Euler beam (adapted from 
[10]). a) The capacitance voltage curve shows two pull-in (A,C) and two release (B,D) instabilities. 
The pull-in voltage at point A corresponds to the transition from no contact (float) to contact with a 
single point (pinned). The pull-in at point C corresponds to the transition from single point contact 
(pinned) to contact over a finite length (flat). The stable equilibria are indicated by solid lines and 
the unstable by dotted lines, b) The displacement shape of the beam at 5 points (A-E) on the C — V 
curve. 



Some finite element packages do provide basic path following methods to solve 
non-hnear problems, however the control over the followed path and the detection 
of stability is less advanced than in the dedicated path following packages ||6l dis- 
cussed in the previous section. To analyze non-linear MEMS systems with path fol- 
lowing methods in a FEM package, we have written a script for Com SOL |3J. The 
script is similar to that presented in Q and predicts the initial condition for the next 
calculation and uses the solver of COMSOL as a corrector The prediction direction 
is a simple linear extrapolation based on the two previous solutions. The predicted 
steplength is based on the previous steplength and the curvature of the path. The 
obtained corrected solution is discarded if its direction or distance is too far from 
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the predicted solution. Since COMSOL provides direct access to the stiffness matrix 
(which is the Hessian of the Hamihonian), the stabiUty of the solutions can be deter- 
mined by checking the sign of the eigenvalues of this matrix. By making use of the 
stiffness matrix, more sophisticated predictor algorithm scripts can be developed. 

In figure|5]an example of a FEM simulation of a voltage actuated single clamped 
beam is shown, which was obtained using the COMSOL script. The results corre- 
spond well to the analytical results in |5i. More details and examples of FEM nu- 
merical path following can be found in | lOj. Although this example is still relatively 
simple, the FEM path following method can theoretically be applied to analyze any 
non-linear system of which the Hamihonian ^ can be written in terms of the ana- 
lytic volume integrals (fTlfTOli. 



9 Conclusions 

It has been shown that path following can be a powerful mathematical tool to de- 
termine the equilibrium solutions of electromechanical systems and their stability 
from the energy and work expressions. The path following of the electromechanical 
equilibrium solutions can be performed by analytical methods, dedicated numeri- 
cal software or finite element method software. Although this paper has focused on 
static systems, the numerical path following technique can also be applied to peri- 
odic dynamic systems (see H), which even extends its usefulness as mathematical 
tool for analyzing electromechanical systems. 
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